November, 1993 



THE QUASI-LINEAR EVOLUTION OF THE DENSITY FIELD 
IN MODELS OF GRAVITATIONAL INSTABILITY 

F. Bernardeau 1 , T. P. Singh 2 , B. Banerjee 2 , S. M. Chitre 2 

1 Canadian Institute for Theoretical Astrophysics 
University of Toronto, Toronto, Ontario, Canada M5S 1A1 
2 Tata Institute of Fundamental Research 
Homi Bhabha Road, Bombay 400 005, India 



ABSTRACT 

Two quasi-linear approximations, the frozen flow approximation (FFA) and the 
frozen potential approximation (FPA), have been proposed recently for studying 
the evolution of a collisionless self-gravitating fluid. In the FFA it is assumed that 
the velocity field remains unchanged from its value obtained from the linear theory 
whereas in FPA the same approximation is made for the gravitational potential. In 
this paper we compare these and the older Zel'dovich approximation by calculating 
the evolution of the density in perturbation theory. In particular we compute the 
skewness, including the smoothing effects, and the kurtosis for the FFA, FPA and 
Zel'dovich approximation and compare their relative accuracy. 

Keywords: galaxies: clustering — cosmology: theory — large-scale structure of 
Universe. 
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1. INTRODUCTION 

Large-scale structures observed today are believed to have developed from small 
density fluctuations generated in the early universe. The growth of these 
fluctuations has been studied by regarding the system to be a collisionless self- 
gravitating fluid. When the amplitudes of the fluctuations are small, perturbation 
theory can be use to study the evolution of the system. In particular the growth rate 
of the rms fluctuations can be described by the linearised equations of fluid motion. 
When the fluctuations eventually become nonlinear, however, the perturbation 
theory is no more accurate and N-body simulations have been widely used to 
overcome this difficulty. The understanding of the physical processes that take 
place in such a self-gravitating fluid necessitates, however, the use of analytical 
models and approximations that can be more easily studied. 

The best known model for describing the mildly nonlinear evolution is due 
to Zel'dovich (1970) (also see Zel'dovich & Shandarin 1989). In this approximation 
the motion of each particle is determined by its initial Lagrangian displacement. A 
good presentation of this approximation has been made by Moutarde et al. 1989 
and by Bouchet et al. 1992 in the frame of a Lagrangian description. The evolution 
of the density field can be studied in this approximation until the formation of 
caustics, when this approach breaks down. Recently two new approximations have 
been proposed, with the aim of improving upon the Zel'dovich approximation - 
(1) the frozen flow approximation (FFA) (Matarrese et al. 1992) and (2) the frozen 
potential approximation (FPA) (Brainerd, Scherrer & Villumsen 1993, Bagla & 
Padmanabhan 1993). In FFA the velocity flows are 'frozen' to their local initial 
linear values, and at any time the velocity of each particle is the one associated to 
the point at which it lies. The evolution of the density is then treated exactly. In 
the second approximation, FPA, the gravitational potential is 'frozen' at its linear 
value. That is, the Eulerian potential is kept constant and the particles obey the 
standard Eulerian equations of motion in this potential. 

All the three approximations mentioned above are naturally consistent with 
the linear theory. However beyond the linear order the density evolution is different 
in each case. To compare them we calculate, assuming Gaussian initial conditions, 
the third and fourth order moments of density by means of perturbation theory. 
We then compare the results obtained for each of these approximations with those 
obtained from perturbation theory using the exact dynamics. 

In Sec. 2 we recall the basic equations of motion and the calculations of 
third and fourth moments of density (Peebles 1980, Fry 1984, Grinstein & Wise 
1987, Bouchet et al. 1992, Bernardeau 1993) for the exact dynamics as well as for 
the Zel'dovich approximation. These computations are repeated for FFA and FPA 
in Sec. 3 and 4 respectively. Throughout we assume an O = 1 spatially flat universe, 
so that the scale factor evolves as a(t) = ao(t/to) 2 ^ 3 . 
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2. BASIC EQUATIONS 

The evolution of a collisionless self-gravitating fluid in an expanding Robertson- 
Walker universe is described by the following equations in the Newtonian limit 
(Peebles 1980), 

s + - v. [(i + <y)v] = o, (i) 

Ct 

v+-v+-[v.V]v=--V0, (2) 
a a a 

V 2 </> = 4nGp b a 2 S. (3) 

Here, 5 = [p(x, t) —p b (t)]/ p b (t) is the enhancement of the true density p(x, t) 
over the mean density Pb(t), v is the proper peculiar velocity, relative to the Hubble 
flow and </>(x, t) is the peculiar gravitational potential. 

From these equations it follows that the density contrast 5(x, t) evolves 
according to the equation 

5 + — 6 - AitGp b 5 = A-KGp b 5 2 + -^-V^.V^ 

a a 2 ^ 

+ ^ ViV,- [(1 + <5)vV] . 

In the linear theory, all the terms on the right hand side of (4) are dropped, 
and 5 = 5^ (x,t) has the solution S^(x,t) = A(x)D(t), with D(t) oc a(t) oc t 2 / 3 
(growing mode) for O = 1. 

It is convenient to define the potential A(x, t) through the relation <fi = 
4irGp b a 2 A : so that V 2 A = 5. The peculiar velocity in the linear theory is 

v« = — VAW, (5) 

where 

A«(x) = --L (M^t. =J^- (6) 
4tt J |x - x'| AnGp b a 2 

The solution for 6 in perturbation theory is obtained via the perturbation expansion 

oo oo oo 

5=^^,v=^vWA = ^AW, (7) 

n=l n=l n=l 
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with respect to the initial Gaussian fluctuations so that 5^ satisfies the equation 
(Fry 1984) 



2d , 



a fc=l 

+ 47rGp 6 V^( fe )V i A (n - fe) + irV l V J v ik)t v^- k)j (8) 

a 2 

1 fc-i 

+ — V v ! vi m) /" m)! t) (n " fc) 



m—l 

The derivation of the behavior of the first cumulants of the density 
distribution at large scale can be done assuming Gaussian initial conditions. It 
can be shown that in general (Goroff et al. 1986, Bernardeau 1992) 

<F> c = S p (S 2 ) p -\ (9) 

where S p is a coefficient that depends weakly on the cosmological parameters. 
Moreover when the smoothing effects are neglected these coefficients are 
independent of the shape of the power spectrum. Bernardeau (1992) gives a method 
to derive the whole series of the coefficients for the exact dynamics, but we consider 
here only the first two coefficients £3 and S4. 

Let us first recall the principle of the calculation of S3 and S4 for the exact 
dynamics. The derivation of the skewness of the distribution function requires the 
calculation of the density contrast at second order, 5^ 2 \ It is obtained by setting 
n = 2 in equation (8), 

g{2) + ^(2) _ 4nGpb6 (2) = 4nGf)b ^(1)2 + Vi(5 <l) V . A (D) 



a 2 



The solution for the growing mode is 

«M = 5«W+tf>A« + HA« a . (10) 



Using this solution, it is straightforward to calculate the skewness £3 = (5 3 )/(5 2 ) , 
which to the lowest order is 3(5 (1 '> 2 5 (2 '>) / (5^ 2 ) 2 = 34/7 (Peebles 1980). 

The calculation of the fourth cumulant involves the knowledge of the density 
field at the third order in perturbative calculation. It can be shown that the 
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coefficient £4 can be written (Fry 1984, Bernardeau 1992) as £4 = 12R a + 414, 
where 

R a = A{5^ 2 5^ 2 )/{5^ 2 )\ (11) 

and 

R b =(6WW>)/(6W)\ (12) 

The value of R a can be easily obtained from Eqn. (10), R a = (34/21) 2 = 2.62. 
The third order term of the density contrast, S^ 3 \ can be calculated more easily in 
Fourier space, by first defining 

*(k) = I J d 3 x5(x)e ik - x , (13) 

and similar transforms for v and A. Eqn. (8) for n = 3 then gives a solution for 
5^ s \ which can be used to show that Rb = 682/189 = 3.61, and that 64 = 45.88 
(see Fry 1984, Bernardeau 1992). 

These results, however, concern the behaviour of the cumulants at a given 
point. When the field is filtered at a given scale, the values of £3 and £4 have 
to be changed. We recall here the results obtained for a top hat window function 
(Juszkiewicz, Bouchet & Colombi 1993, Bernardeau 1993). They read, 

34 , . 

S 3 = y+7i (14) 

and 

. 60712 62 7 2 2 

* 4 = W + T 7l+ 3 7l + 3 72 (15) 
where 71 and 72 are the first two logarithmic derivatives of the variance with scale, 

7 1 = ^1 p ( 16 ' 

d log R 

dlog R 

The derivation of these smoothing effects is based on geometrical properties of the 
top hat window function given by Bernardeau (1993). 

For the Zel'dovich approximation the behaviour of the cumulants at large 
scale is similar to the one encountered in the real dynamics but the coefficients 
S p are slightly changed due to the approximation that is made (Grinstein & Wise 
1987). Recently Bernardeau (1993) has derived the expression of the coefficients S3 
and £4 when a top hat filter is applied to the density field, 

Sl el = A + lu (18) 
nZel 272 50 7 2 2 

St 1 = ^- + Y ll + 3 71 + 3 72 ' (19) 
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This quantitative change of the large-scale cumulants is general for any 
approximative dynamics starting with Gaussian initial conditions. The coefficients 
S p then turn out to be a good tool to test the various approximations with each other 
by comparing the values of these coefficients. In the next parts we will generalize 
the results obtained for the Zel'dovich approximation to the frozen flow and the 
frozen potential approximations. 



3. FROZEN FLOW APPROXIMATION 

The frozen flow approximation (FFA) which was proposed by Matarrese et al. (1992) 
is best defined using slightly different variables from those in Eqns. (1-3). Using the 
scale factor 'a' as the time variable, define the comoving peculiar velocity u = dx/da 
= v/ad. 77 = 1 + 5, ip = (3tg/2ag)0, where a(t) = a (t/t ) 2/3 . Eqns. (1-3) then 
reduce to 

^ +?? V.u = 0, (20) 
da 

du 3 3 „ , . 

da 2a la 

V 2 ^ = S/a, (22) 

where d/da = d /da + u.V. 

FFA is defined by assuming that the velocity field u is steady: du/da = 0; 
that is, stream lines are frozen to their initial shape. The frozen value of u would 
then be the constant value it has in the linear theory: 

uffa(x) = -VV'ljjv(x). (23) 

(This of course implies that v = vlin = v^, as given by Eqn. (5).) The uffa of 
Eqn. (23) is a solution of the Euler equation (21) provided ip is approximated to be 

CI 

•0FF^(x, t) = IpLIN - g (V^ L /jv) 2 . (24) 
In the notation of Eqns. (1-3), FFA corresponds to 



3 \2aQ 



4 ^...2/„ A (i)V (25) 



0« - ^Gp b a* (VAW) 
4>M + 0( 



The form of (p^ implies that it is second order in perturbation, but it is not the 
same as the second order potential in true non-linear evolution. To study the density 
evolution in FFA, we first note that Eqns. (1-3) give 

6 + — 6 = \(1 + 5)V 2 <P + -1V(5.V0 

a a z a z ^6) 

+ ^ ViWj [(1 + 8)v l v j ] . 

In here, we substitute for v as vffa and for <p as 4>ffa- Next, we implement the 
perturbation expansion of Eqn. (7) to get the following equation for #( 2 ): 



a a z a A ^7) 

a 1 L J 

This equation has the solution 

»f FA = \^ + (28) 

which should be compared with the 5^ for the true evolution, in Eqn. (10). 
From here, it is straightforward to carry through the calculation of skewness as 
in Peebles (1980), Section 18, since the only change in the solution 5^ is that in 
the coefficients. The result is S^ FA = 3, as compared to the true value of 34/7. 
The smoothing effects on £3 can be easily calculated and the final result for a top 
hat window function reads, 

S[ FA = S+lL. (29) 

To obtain £4 in FFA, we first find R a from Eqn. (11) by substituting the 

(2) 

solution S FFA . This gives, upon angle averaging, as in Fry (1984), R a = 1.0. From 
(26), the Eqn. for 5^ in FFA reads 

#3) + 2o#3) = 1 r 5 (2) v 2 (l) +5 (1) V 2 (2)| 

+ JL {v^.V^ + V^.V^} (30) 

a 2 L J 
Before doing the Fourier transform, we note that 

V^ 2 ) = -^-Gp b a 2 (vA«.v) VA«, (31) 
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VV = -|Gp t a 2 {A<»^» + A»»A«}. (32) 
Using the transform (13) and its inverse 

*(x) = J^J d s k~5(k)e^ : (33) 

the Fourier transform of a product can be written as a convolution (Fry, 1984): 



FT{F 1 (x)---F N (x)} 



N 



i r d 3 h d 3 k 

W J (2^)3 "' (2^)3 



F 1 (ki)---Fjv(k A r) = Fi* 
When applied to Eqn. (30) for 8pp A , this gives 



FM^ki-k) 

• • * F N . 



(34) 




a 2 

The solution can be written down in analogy with Eqn. (46) of Fry (1984), so we 
do not put it down explicitly, except to note that now 5pp A of Eqn. (28) should 
be used. Using this 5^ 3 ^ and carrying out angular averaging as in Fry's paper gives 
Rb = 1, and hence that 

S[ FA = 16. (36) 

The derivation of the smoothing effects for S4 is slightly more complicated and will 
not be given. 
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4. FROZEN POTENTIAL APPROXIMATION 

The frozen potential approximation (FPA) was proposed by Brainerd et al. (1993) 
and by Bagla & Padmanabhan (1993). FPA is defined by keeping the potential <j) 
constant at its value ^^(x) in the linear theory, so that (p^ and higher terms in 
the perturbation expansion of (p are set to zero. However, unlike in FFA, v is not 
approximated, but is to be obtained from the Euler equation (2) with </>(x, t) = <j)d> 
(x). Thus Eqn. (26) for the evolution of density, when written for the case of FPA, 
becomes 

s + = -1(1 + a)vy ^ + i^vs.v^ 

a or or 

+ ^ViVj [(i + o>v] . 

Using the perturbation expansion (7), the equation for 5^ is found to be 

#2) + 24^(2) = AnG Pb 5^ 2 + Iv^^.V^ 1 ) 

a 1 a (38) 



+ — ViVj 
a 2 



The solution is 

4 P) = + 1^>A« + 1a!» 2 , (39) 

as contrasted to the true #( 2 ) in Eqn. (10), and Sp FA in Eqn. (28). Once 
again, the skewness is easy to work out, following Peebles (1980), and the result is 
S[ PA = 17/5. The derivation of the smoothing effects gives 

= % + ^7, (40) 



The calculation of the kurtosis of the density field can be done as usual: 
Substituting this expression for o^ 2 ) in (11) and applying angle averaging as before, 
gives R a = (17/15) 2 = 1.28. The equation for 5^ in FPA is 



#3) + ^(3) = 4nGpb6 (l) 6 (2) + J_ V5 (2)_ V0 (1) 

a or 



S (%Wi v M3 +v Wi v V)3 +v Wi v Wi 



(41) 



as contrasted to the equation (30) for Sp FA . Here, 5^ is the FPA solution Eqn. 
(39), and the solution for v^ 2 ^ can be found from the continuity equation to be 



v -_1^VA (2) -tf^M 1 ) 

W FPA — FPA ' 



r( 2 ) 



(42) 
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Eqn. (41) is a special case of the equation for the true S^ 3 \ and has been obtained 
simply by setting tp^ = 0. It is thus easier to handle than the equation for Sp FA . 
Carrying out the Fourier transform and angular averaging precisely as in Fry (1984) 
gives Rb — 457/315 = 1.45, and hence 

S( PA = 21.22. (43) 



5. CONCLUSIONS 

In Table 1, we display the leading order mean values of the third and fourth moments 
for the true and approximate evolution. For the third moment, the smoothing 
corrections for a top hat window function are included, (n+3 = — d log (5 2 )/d log R) . 

Table 1 



Third moment, £3 Fourth moment, £4 



True evolution 4.86 - (n + 3) 45.88 

Zel'dovich 4.00 - {n + 3) 30.22 

Frozen potential 3.40 - 0.7 (n + 3) 21.22 

Frozen flow 3.00 - 0.5 (n + 3) 16.00 



How do our results compare with those from simulations using N-body, FFA 
and FPA? (For results on these simulations, see Matarrese et al. 1992, Brainerd 
et al. 1993, Bagla & Padmanabhan 1993, Mellott et al. 1993). At first sight, one 
might think that since both FFA and FPA simulations are more similar to N- 
body than Zel'dovich approximation, our analytical results conflict with results of 
density evolution from the simulations. However, while FFA and FPA both prevent 
the thickening of pancakes, they do not do well in moving mass to the right place. 
Brainerd et al. (1993) find from their analysis that the cross-correlation between 
N-body and Zel'dovich is higher than that between N-body and FPA. Also, Mellott 
et al. (1993) find that FFA does poorly in cross-correlation with N-body. 

At this stage it is useful to compare the second order solution for the density 
field in the various approximations. For the Zel'dovich case, the second order 
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solution is (Bouchet et al. 1992) 

jM^iW+tf'AW + (44) 

and the second order solution for FFA and FPA is given in Eqns. (28) and 
(39) respectively. The same terms appear in the three approximations, but the 
coefficients are different. The coefficients in the Zel'dovich case are the closest to 
the true density evolution, followed by FPA and then FFA, and this feature is 
clearly reflected in the respective values of the skewness. It is now interesting to 
try to understand why these approximations underestimate such a quantity. The 
skewness can be seen as a measure of the ability of the system to create rare dense 
spots. This interpretation seems to be confirmed by the relation of the spherical 
collapse with the value of the skewness (see Bernardeau 1993). As can be seen in 
Matarrese et al. (1992), the FFA is less accurate than the Zel'dovich approximation 
(and the real dynamics) for the spherical collapse: the acceleration is too weak to 
concentrate the matter efficiently. This result is again borne out by the result on 
the kurtosis. 

When the smoothing effects are taken into account, however, the results for 
the skewness seem to attenuate these effects. In particular for n = — 1 the Zel'dovich 
approximation, the FFA and the FPA all give the same result, S3 = 2 (instead of 
the exact value £3 = 2.86). Actually the smoothing corrections are sensitive to the 
tidal effects in the density field (they come from the term 6 pA^) and only the 
Zel'dovich approximation gives the right coefficient for this term. The other two 
approximations underestimate the tidal effects. The FFA even fails to give a term 
containing a quadrupole contribution (in A^A^ 1 - ). This is a major consideration 
for a practical use of these approximations. The disruption of objects, for instance, 
is expected to be less accurate in the FFA or the FPA than in the real dynamics. 

FFA was proposed to improve upon the Zel'dovich approximation, by 
avoiding shell-crossing. In a similar spirit, FPA attempts to improve upon FFA 
and Zel'dovich, by keeping only the potential linear, but evolving both density 
and velocity exactly (shell-crossing does take place). However, our results suggest 
that FFA and FPA are only partially successful in their aim, and the quasilinear 
regime is rather poorly approximated by these approximations. Thus while 
the simulations and our analytical results are two different ways of testing the 
various approximations, both the means indicate the need for a more careful 
comparison, before FFA and FPA can be adopted as improvements over the 
Zel'dovich approximation. These approximations might of course be interesting 
for analytical studies, but one should generally exercise caution in their use. 

F.B. and T.P.S. would like to thank the organizers of the School on Cosmology and 
Large Scale Structure, at Les Houches, where some of these ideas were discussed. 
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